Analytic study of clustering in shaken granular material using zero-range processes 
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We show that models used to described granular clustering due to vertical shaking belong to the 
class of zero-range processes. This correspondence allows us to derive exactly in a very easy and 
straightforward manner a number of properties of the models like particle distribution functions, 
phase diagram, and characteristic time of clusterization. 
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I. INTRODUCTION 

Non-equilibrium phase transitions were observed in 
many simple systems |il LI |3f . Recently, it was also found 
that shaken granular material exhibits clustering depend- 
ing on the shaking strength Q . In his paper Jens Eggers 
suggested a model for the description of the clustering of 
vertically shaken granular material. Originally, it was in- 
troduced for the two box setup. Later it and its modified 
version were used to describe experiments having more 

boxes H US SUES EH 

The analytic studies of the above mentioned papers 
were difficult and in many cases did not give general re- 
sults. In the case of exclusion models, it was proved 
E3 that the correspondence to an already solved model, 
namely the zero-range process [Tl| can be of great help, 
as many results can be obtained directly. 

The aim of this paper to derive analytically for the 
general case the probability distribution, phase diagram 
and the coarsening time using the above correspondence. 
We note that specific cases with small number of boxes 
are already solved H H H U 03 G3 



II. THE MODEL 



depend only on the number of particles on the site of 
departure and external parameters but independent of 
the properties of the target site. Here, we consider only 
the symmetric case when the hopping to the left or to 
the right is equally probable. 

The important attribute of the zero-range process is 
that it yields a steady state described by a product mea- 
sure. By this it is meant that the steady state prob- 
ability P({n^}) of finding the system in configuration 
{nx,n2, ■ ■ -Ul} is given by a product of factors 
that are called marginals 
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where Z(L, N) is the normalization factor. For the zero- 
range process fin^) is given by 
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The experimental setup looks as follows. The sys- 
tem consists of a container separated into boxes by walls 
which are open upwards. The whole system is then ver- 
tically shaken and particles can hop above the walls from 
one compartment to the other. The model of Eggers |j| 
defines particle fluxes between the boxes based upon the 
physical properties of the system like shaking strength 
and local particle density, etc.. Two different steady 
states were found, both in the experiments and in the 
model: a homogeneous where the containers hold roughly 
equal number of particles, and a condensed steady state 
where one compartment contains nearly all particles. 

We first define the zero-range process following E3 : Q 
We consider a one-dimensional finite lattice of L sites and 
periodic boundary conditions. The total particle number 
is denoted by N. 

The dynamics of the system is given by rates u(n), at 
which a particle leaves a site. The hopping rates u(n) 



The marginals are defined up to a multiplicative factor. 

It is important to note that the steady state particle 
distribution can be formally read off generally for all N 
and L. This will allow us to calculate the general phase 
diagrams. 

The Eggers model !j| has a quadratic temperature de- 
pendence while Lipowski and Droz y| uses a linear form. 
For all quantities we use an index E for Eggers model and 
LD for the Lipowski-Droz model. The fluxes are defined 
in the following way 
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exp(-l/[T + A(l-n/JV)]), (3) 



where A, B, To, and A contain physical quantities. It 
is of importance that B is proportional to the inverse 
shaking strength while A increases with the strength of 
the agitation. The u(n) functions are compared on Fig. 



Our main purpose in this paper is to study the clus- 
terization process of this model. Furthermore, as cases 

with L = 2,3HH|l|l|!E3,E3are already SOlVed W6 
focus on the general case with at least moderate L and 

large N 3> 1 particle number. 
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FIG. 1: Sample hopping probabilities. Solid line E dashed 
LD model. 



IV. PHASE DIAGRAM 

We already discussed that both models have two pos- 
sible stable steady states which give rise to four different 
cases on the phase diagram: I) only the symmetric state 
is stable, II) only the asymmetric state is stable, III) both 
states are stable but asymmetric is more probable, IV) 
both states are stable but symmetric is more probable. 
In this section we determine the three lines separating 
these states. 

On line III-IV the symmetric solution loses its stabil- 
ity. In the symmetric steady state all sites contain N/L 
particles. Only a flux curve with positive derivate may 
be stable for many sites thus N/L must be less than the 
position of the maximum of the flux curve. This gives 
the following two relations 



III. PARTICLE DISTRIBUTION 



(E) B = L 2 /N 2 
(LD) T = Ja/L- 



A(l-1/L). 



(6) 



Before calculating the general particle distribution 
function we remark that in both models only two dif- 
ferent steady states may be stable: (i) The symmetric, 
where all sites contain the same amount of particles, and 
(ii) The asymmetric, or condensed steady state, where 
one site contains more particles than the others which 
have equal number of particles. 

The reason why only those two are stable is that in 
the steady state all sites must have the same flux. The 
hopping probability curve has a maximum thus sites may 
be on both sides of the maximum at the positions with 
equal flux. However, when the flux has negative derivate 
only one site may be at a given position because more 
would be unstable against fluctuations. 

Thus we calculate only the probability distribution 
function for the case where L — 1 sites contain the 
same number of particles a the remaining site thus holds 
b = N — (L — l)a particles. The marginals take the form 
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We are lucky that for both models the unnormalized 
probability distribution can be obtained in a closed form: 



P E ({a,...a,b}) = 
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The other two line can only be determined through the 
probability distribution function. 



A. The LD model 

Let us define the particle ratios: x = a/N, f3 = b/N . 
The following equation describes the extremal points of 
the probability distribution function 
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A(l-/3) T +A(l-x)' 



(7) 



It is obvious to see that (3 — x is always a solution and 
that the above equation is N independent. The exact 
position of the clusterized state can only be determined 
numerically. 

An important question is whether the transition from 
symmetric to asymmetric steady state is continuous or 
not. The third order Taylor expansion of Eq. J7J) around 
the symmetric solution of a; = 1/L shows that if L > 2 no 
continuous transition is possible as there is no solution 
of the equation. While for L = 2 continuous solution is 
possible if A < 2/3. This result was also shown with 
a general argument in ^(j- It is important to reiterate 
because in the following calculations are valid for x -C I 
which can only be used in case of discontinuous transi- 
tions. 

The line I-IV, the limit of the stability of the asym- 
metric steady state can be derived from Eq. Q. As- 
suming that i<1 and the fact that at the point where 
other roots except from [3 — x appear the value and the 
derivate of both sides of the equation should be equal. 
This results in the following implicit equation 
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FIG. 2: Phase diagrams. Top (a) LD for L = 6 bottom (b) 
model E. 



The line III-IV, can not be exactly determined because 
we do not know exactly the roots of Eq. Q . Because of 
the previous arguments a good estimate can however, be 
found by assuming x — and (3=1 which is a very good 
approximation in most parts of the phase space. The 
relation P L d(0, ... 0, TV) = P L d{N/L, . . . N/L) gives an 
implicit equation for A and Tq 
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An example phase diagram for L = 6 can be seen on 
Fig. a. 



B. The E model 

Following the same route as for the LD model we derive 
the equation for the position of the steady states in the 
leading order of N. 



(P 2 



)BN 



2ln((3/x) + 0(l/N). (10) 



The N dependence of the equation disappears if one 
chooses B oc N~ 2 . It is also obvious that the above 
equation has the solution of x — (3. Other solutions must 
satisfy x < [2(L — 1)]~ 1 . This does not allow a continu- 
ous phase transition for L > 2 but continuous transition 
always exists for L — 2 as was shown in |5j- 



Line I-IV can be determined from the above results 

BN 2 = 81nL + 0(1). (11) 

Line III-IV can be determined from the equality of the 
probability distribution functions as 



BN 2 =61nL + 0(l). 



(12) 



The above two equations suffer from strong corrections 
to the leading order InL term. However, this accuracy 
is enough the get a clear picture of the phase diagram 
which is shown on Fig. 0b. 



V. COARSENING TIME 

The coarsening process in general drives the system to 
have bigger but a smaller number of macroscopic clusters 
until only one prevails. As the clusters become bigger 
the coarsening process slows down enormously because 
the size of the clusters to be dissolved is bigger and the 
relative flux is much smaller. Thus, in general it is enough 
to study the competition of the last two big clusters [Tj| 
to get the asymptotic behaviour of the coarsening time 

T. 

The above picture applies only if the if the homoge- 
neous steady state is not stable. In the opposite case, 
where the symmetric steady state is also stable a fluctu- 
ation large enough is needed to drive the system out of 
the potential well for which the time needed is in general 
proportional to exp(— A^ 2 ). Here we study only the case 
when the system is in section II of the phase space. 

For obvious reasons, the two big clusters are generally 
at ~ L/2 distance. If a big cluster looses a particle to 
its neighbour the particle starts a random walk. If the 
distance between the two clusters is L/2 then the prob- 
ability of reaching the other is 2/L and it takes T 2 /4 
timesteps. Thus the coarsening time can be written as 



N/2-xi 
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u{N/2-n)-u(N/2 + n) 



(13) 



The coarsening time is be dominated by the state when 
both clusters have almost the same number of particles. 
Thus, we use the first order Taylor series around this 
point: 
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For the E model we have to identify the N dependence 
of A. This parameter sets the speed of the hopping but 
in all previous calculations it was unimportant. We show 
that in the transition regime it should scale as 1 /N 2 ex- 
actly as B. One reason could be the physical derivation 
of Eggers , where one can look up the correspondence 
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FIG. 3: (a) Top: Coarsening time vs. particle number. Box 
number L = 2, 3, 5, 10, 20 for both models scaled together. 
The solid line is iVln(AT/2). The stronger deviations for small 
particle numbers are for the E model. 

(b) Bottom: Coarsening time vs. box number. System sizes 
N = 80, 160, 320 are scaled together for both models. Param- 
eters were chosen as B oc L 2 , To oc d oc 1/L, The solid line is 
a linear. Curves with higher deviations for small box numbers 
are for the LD model and lnr is plotted. For the E model 
\/ln td is shown. 

of A and B, and find that they have the same N de- 
pendence. The other reason is that we want to have the 
same order of fluxes in different systems. The value of 
the maximum of u(n) is AN 2 so it is obvious to have 



A oc l/N 2 in the transition regime. Finally, we can get 
the same N scaling as for the LD model namely 

t oc 2Vm(AT/2), (15) 

which is shown on Fig. U3 

The L dependence of the coarsening time can also be 
read off. For the E model in region II of the phase space 
B has a L dependence faster than 1? . It is easy to see 
that if B oc L 1 then the coarsening time grows as te oc 
exp(T 7 ). For the LD model the dominating term will 
also be the exponential as in region II of the phase space 
both A and To should decrease at least as 1 /L for which 
tld oc exp(L). The results are tested against numerical 
simulations on Fig. b. 

As we already mentioned at the model definition B is 
inversely proportional to the shaking strength, while A 
is proportional to the agitation. From physical point of 
view the L scaling of the coarsening time is also the same 
for both models in spite of all differences. 

VI. CONCLUSION 

In this paper we showed that models proposed to de- 
scribed clustering of shaken granular material belongs to 
the group of zero-range processes which is very efficient 
in describing non-equilibrium phase transitions. 

The probability distribution which describes the 
steady state of the system can be read off in a straight- 
forward manner. This allows us to determine the phase 
diagram analytically for both models in the general case, 
for L,N < 1. 

The coarsening time is also calculated for both models 
and tested against numerical simulations. The two mod- 
els have qualitatively different phase diagram but sur- 
prisingly similar coarsening time scaling. The particle 
number scaling was found to be the same linear, the box 
number dependence was found to be exponential. 

Part of this work was supported by OTKA F047259 
and by the German-Hungarian Cooperation Fund. 
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